All runs Ulva & Hypnea rETRmax Analysis, Script Chunks, and Plots

This is the analysis of all the Ulva lactuca salinity and nutrient experiments conducted on the lanai in St. John 616 from September 2021 to October 2022. These experiments incorporated four paired salinity and nutrient treatments with three temperatures. Each of the first four runs produced an n = 2 and was repeated initially 8 times for a total of n = 16. Data gaps were identified and filled in February, April, and October 2022. This output reflects all data totaling five treatments for Ulva lactuca and six for Hypnea musciformis.

Packages loaded:

library(lme4)
library(lmerTest)
library(afex)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load this dataset normalized to quantum efficiency of photosynthesis per Silsbe and Kromkamp 2012

all_runs_photosyn_data <- read.csv("../data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")

Assign several variables as factors

all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)
all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)
all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))
all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)

Subset the species for output. Use Day 9 for final analysis. Remove treatment 2.5 from Ulva dataset

ulva <- subset(all_runs_photosyn_data, Species == "ul" & RLC.Day == 9 & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

#ULVA rETRmax________________________________________________________________

Make a histogram of the data

ulva %>% ggplot(aes(rETRmax)) +
  geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Run model without interaction between the treatments and temperature

ulva_retrmax_model <- lmer(formula = rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)

Make residual plots of the data for ulva

hist(resid(ulva_retrmax_model))

plot(resid(ulva_retrmax_model) ~ fitted(ulva_retrmax_model))

qqnorm(resid(ulva_retrmax_model))
qqline(resid(ulva_retrmax_model))

Check the performance of the model, make a table, plot effects, and get r2

performance ::check_model(ulva_retrmax_model)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(ulva_retrmax_model)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.5925302 0.6884433
summary(ulva_retrmax_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 1826.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8395 -0.6155  0.0792  0.5149  3.4121 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  21.987   4.689  
##  Run       (Intercept)   6.911   2.629  
##  RLC.Order (Intercept)   5.089   2.256  
##  Residual              110.401  10.507  
## Number of obs: 240, groups:  Plant.ID, 95; Run, 7; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)     34.821      2.989  9.287  11.648 7.53e-07 ***
## Treatment1      20.947      3.229  7.232   6.486 0.000295 ***
## Treatment2      22.032      3.229  7.232   6.822 0.000214 ***
## Treatment3      38.154      3.229  7.232  11.815 5.47e-06 ***
## Treatment4      39.478      3.229  7.232  12.225 4.32e-06 ***
## Temperature27   -4.440      2.497 30.608  -1.778 0.085323 .  
## Temperature30   -2.153      2.285 66.092  -0.942 0.349434    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.644                                   
## Treatment2  -0.644  0.779                            
## Treatment3  -0.644  0.779  0.779                     
## Treatment4  -0.644  0.779  0.779  0.779              
## Temperatr27 -0.399  0.003  0.003  0.003  0.003       
## Temperatr30 -0.387  0.000  0.000  0.000  0.000  0.475
tab_model(ulva_retrmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  rETRmax
Predictors Estimates std. Error CI Statistic p df
(Intercept) 34.82 2.99 28.93 – 40.71 11.65 <0.001 229.00
Treatment [1] 20.95 3.23 14.58 – 27.31 6.49 <0.001 229.00
Treatment [2] 22.03 3.23 15.67 – 28.40 6.82 <0.001 229.00
Treatment [3] 38.15 3.23 31.79 – 44.52 11.81 <0.001 229.00
Treatment [4] 39.48 3.23 33.11 – 45.84 12.22 <0.001 229.00
Temperature [27] -4.44 2.50 -9.36 – 0.48 -1.78 0.077 229.00
Temperature [30] -2.15 2.28 -6.65 – 2.35 -0.94 0.347 229.00
Random Effects
σ2 110.40
τ00 Plant.ID 21.99
τ00 Run 6.91
τ00 RLC.Order 5.09
ICC 0.24
N Run 7
N Plant.ID 95
N RLC.Order 6
Observations 240
Marginal R2 / Conditional R2 0.593 / 0.688 
plot(allEffects(ulva_retrmax_model))

Construct null model to perform likelihood ratio test REML must be FALSE

ulva_retrmax_treatment_null <- lmer(formula = rETRmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_retrmax_model2 <- lmer(formula = rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_retrmax_treatment_null, ulva_retrmax_model2)
## Data: ulva
## Models:
## ulva_retrmax_treatment_null: rETRmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_retrmax_model2: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                             npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_retrmax_treatment_null    7 1983.9 2008.2 -984.94   1969.9          
## ulva_retrmax_model2           11 1871.1 1909.4 -924.57   1849.1 120.74  4
##                             Pr(>Chisq)    
## ulva_retrmax_treatment_null               
## ulva_retrmax_model2          < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_retrmax_temperature_null <- lmer(formula = rETRmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_retrmax_model3 <- lmer(formula = rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_retrmax_temperature_null, ulva_retrmax_model3)
## Data: ulva
## Models:
## ulva_retrmax_temperature_null: rETRmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_retrmax_model3: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                               npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_retrmax_temperature_null    9 1870.0 1901.3 -925.99   1852.0          
## ulva_retrmax_model3             11 1871.1 1909.4 -924.57   1849.1 2.8537  2
##                               Pr(>Chisq)
## ulva_retrmax_temperature_null           
## ulva_retrmax_model3               0.2401

Plots

ulva %>% ggplot(aes(treatment_graph, rETRmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "rETRmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 170) + stat_mean() + 
  scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
  geom_hline(yintercept=0, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Summarize the means for rETRmax

ulva %>% group_by(Treatment) %>% summarise_at(vars(rETRmax), list(mean = mean))
## # A tibble: 5 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          32.6
## 2 1          53.2
## 3 2          54.3
## 4 3          70.4
## 5 4          71.8
ulva %>% group_by(Run) %>% summarise_at(vars(rETRmax), list(mean = mean))
## # A tibble: 7 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      64.7
## 2 2      58.4
## 3 3      64.6
## 4 3.5    66.3
## 5 4      61.2
## 6 6      35.8
## 7 6.5    29.4
ulva %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(rETRmax), list(mean = mean))
## # A tibble: 5 × 3
## # Groups:   Treatment [5]
##   Treatment RLC.Day  mean
##   <fct>       <int> <dbl>
## 1 0               9  32.6
## 2 1               9  53.2
## 3 2               9  54.3
## 4 3               9  70.4
## 5 4               9  71.8

Add growth rate to the dataset and subset by species

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$lunar.phase <- as.factor(growth_rate$lunar.phase)

Make a new column for weight change (difference final from initial)

growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
ulva$lunar.phase <- (gr_ulva$lunar.phase)

Plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_etr_graph <- ggplot(ulva, aes(x=rETRmax, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca rETRmax vs Growth Rate", x = "rETRmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph
## `geom_smooth()` using formula = 'y ~ x'

HYPNEA____________________________________________________________________________________________________________

There was no D9 RLC for hm6-4 on 11/12/21 but had to remove hm6-4 from 10/9/21 below to match growth data because replicate appeared dead

hypnea <- subset(all_runs_photosyn_data, Species == "hm" & RLC.Day == 9 & uid != "2021-10-09_hm6-4")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

For Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) and hm6-4 on 10/9/21 because it was white and also white and dead hm1-1 on 10/12/22 and hm1-2 on 4/29/22 causing issues of influential observations hm1-2 for both rETRmax and Ek – leaving them in dataset because no good reason to believe not good data

gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

Make a histogram for hypnea

hist(hypnea$rETRmax, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)

hypnea %>% ggplot(aes(rETRmax)) +
  geom_histogram(binwidth=5, fill = "#7D0033", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

Run model for rETRmax

hyp_retrmax_model <- lmer(formula = rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea)

Make residual plots of the data

hist(resid(hyp_retrmax_model))

plot(resid(hyp_retrmax_model) ~ fitted(hyp_retrmax_model))

qqnorm(resid(hyp_retrmax_model))
qqline(resid(hyp_retrmax_model))

Check the performance of the model, get r2, make table, and plot effects

performance::check_model(hyp_retrmax_model)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(hyp_retrmax_model)
##            R2m       R2c
## [1,] 0.2758676 0.6620738
summary(hyp_retrmax_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: hypnea
## 
## REML criterion at convergence: 2279.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1653 -0.4931 -0.0716  0.3862  4.0746 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  97.002   9.849  
##  Run       (Intercept)  35.742   5.978  
##  RLC.Order (Intercept)   4.414   2.101  
##  Residual              120.012  10.955  
## Number of obs: 286, groups:  Plant.ID, 143; Run, 8; RLC.Order, 6
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    67.24266    5.12098   5.80941  13.131 1.54e-05 ***
## Treatment1     -0.04109    5.85801   5.16221  -0.007    0.995    
## Treatment2     -2.01297    5.85801   5.16221  -0.344    0.745    
## Treatment2.5    9.99208    7.91551   3.76075   1.262    0.279    
## Treatment3      3.72557    5.85801   5.16221   0.636    0.552    
## Treatment4      1.75044    5.86887   5.20152   0.298    0.777    
## Temperature27 -18.99828    3.13838  11.86996  -6.054 6.00e-05 ***
## Temperature30 -19.45096    2.98374  24.17631  -6.519 9.31e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.747                                          
## Treatment2  -0.747  0.927                                   
## Treatmnt2.5 -0.552  0.483  0.483                            
## Treatment3  -0.747  0.927  0.927  0.483                     
## Treatment4  -0.746  0.925  0.925  0.482  0.925              
## Temperatr27 -0.299  0.003  0.003  0.000  0.003  0.003       
## Temperatr30 -0.294  0.000  0.000  0.000  0.000  0.003  0.487
tab_model(hyp_retrmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  rETRmax
Predictors Estimates std. Error CI Statistic p df
(Intercept) 67.24 5.12 57.16 – 77.32 13.13 <0.001 274.00
Treatment [1] -0.04 5.86 -11.57 – 11.49 -0.01 0.994 274.00
Treatment [2] -2.01 5.86 -13.55 – 9.52 -0.34 0.731 274.00
Treatment [2.5] 9.99 7.92 -5.59 – 25.58 1.26 0.208 274.00
Treatment [3] 3.73 5.86 -7.81 – 15.26 0.64 0.525 274.00
Treatment [4] 1.75 5.87 -9.80 – 13.30 0.30 0.766 274.00
Temperature [27] -19.00 3.14 -25.18 – -12.82 -6.05 <0.001 274.00
Temperature [30] -19.45 2.98 -25.32 – -13.58 -6.52 <0.001 274.00
Random Effects
σ2 120.01
τ00 Plant.ID 97.00
τ00 Run 35.74
τ00 RLC.Order 4.41
ICC 0.53
N Run 8
N Plant.ID 143
N RLC.Order 6
Observations 286
Marginal R2 / Conditional R2 0.276 / 0.662 
plot(allEffects(hyp_retrmax_model))

Construct null model to perform likelihood ratio test REML must be FALSE

hypnea_retrmax_treatment_null <- lmer(formula = rETRmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_retrmax_model2 <- lmer(formula = rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_retrmax_treatment_null, hypnea_retrmax_model2)
## Data: hypnea
## Models:
## hypnea_retrmax_treatment_null: rETRmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_retrmax_model2: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                               npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_retrmax_treatment_null    7 2334.6 2360.2 -1160.3   2320.6          
## hypnea_retrmax_model2           12 2334.8 2378.7 -1155.4   2310.8 9.8341  5
##                               Pr(>Chisq)  
## hypnea_retrmax_treatment_null             
## hypnea_retrmax_model2            0.08007 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hypnea_retrmax_temperature_null <- lmer(formula = rETRmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_retrmax_model3 <- lmer(formula = rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_retrmax_temperature_null, hypnea_retrmax_model3)
## Data: hypnea
## Models:
## hypnea_retrmax_temperature_null: rETRmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_retrmax_model3: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                                 npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_retrmax_temperature_null   10 2363.5 2400.1 -1171.8   2343.5          
## hypnea_retrmax_model3             12 2334.8 2378.7 -1155.4   2310.8 32.742  2
##                                 Pr(>Chisq)    
## hypnea_retrmax_temperature_null               
## hypnea_retrmax_model3            7.765e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plots

hypnea %>% ggplot(aes(treatment_graph, rETRmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="Treatment", y= "rETRmax (μmols electrons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 170) + stat_mean() + 
  scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
  geom_hline(yintercept=0, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Plot a regression between the photosynthetic independent variables of interest and growth rate

hypnea_growth_etr_graph <- ggplot(hypnea, aes(x=rETRmax, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis rETRmax vs Growth Rate", x = "rETRmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor(label.y = 150)
hypnea_growth_etr_graph
## `geom_smooth()` using formula = 'y ~ x'

Summarize the means

hypnea %>% group_by(Treatment) %>% summarise_at(vars(rETRmax), list(mean = mean))
## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          54.4
## 2 1          55.3
## 3 2          53.3
## 4 2.5        64.4
## 5 3          59.0
## 6 4          57.3
hypnea %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(rETRmax), list(mean = mean))
## # A tibble: 6 × 3
## # Groups:   Treatment [6]
##   Treatment RLC.Day  mean
##   <fct>       <int> <dbl>
## 1 0               9  54.4
## 2 1               9  55.3
## 3 2               9  53.3
## 4 2.5             9  64.4
## 5 3               9  59.0
## 6 4               9  57.3